function eq_mit = solve_eq_mit(sim, eq_SS, param, glob, options)

    diff = 10;
    %iter = 0;
    
    T_plot = 20;
    
    while (diff > 1e-8)
        
        sim_out  = solve_bf(sim, eq_SS, param, glob, options);
        
        diffC = max(max(abs(sim_out.C - sim.C)/abs(sim.C)));
        diffW = max(max(abs(sim_out.W - sim.W)/abs(sim.W)));
        
        diff  = max([diffC diffW])
        
        Labor       =  sim_out.l;
        total_labor =  sim_out.L; sum(sim_out.mu.*Labor);
        
        figure(47)
        
        subplot(3,3,1)
        hold off
        plot([sim_out.C - 1]); drawnow; title('C'); drawnow;
        xlim([1 T_plot]);
        
        figure(47)

        subplot(3,3,2)
        hold off
        plot(sim_out.D - eq_SS.D); drawnow; title('D'); drawnow;
        xlim([1 T_plot]);

        figure(47)
        
        subplot(3,3,3)
        hold off
        plot([total_labor]); title('Total labor demand'); drawnow;
        xlim([1 T_plot]);

        figure(47)
        
        mass = sum(sim_out.mu);
        
        subplot(3,3,4)
        hold off
        plot(mass/mass(end)); title('Mass of firms'); drawnow;
        xlim([1 T_plot]);
        
        figure(47)
        subplot(3,3,5)
        hold off
        plot(sum(sim_out.exit.*sim_out.mu)./sum(sim_out.mu)); title('Exit rate'); drawnow;
        xlim([1 T_plot]);

        figure(47)
        subplot(3,3,6)
        hold off
        plot(sim_out.ve); title('Value of entry'); drawnow;
        xlim([1 T_plot]);

        costs = sim_out.l;
        agg_prods = repmat(sim.A', glob.Nsf,1);
        prods = repmat(exp(glob.sf(:,2)), 1, options.T);
        prods = agg_prods.*prods;
        wages = reshape(sim_out.W,1,[]);
        wages = repmat(wages, glob.Nsf, 1);
        mkps  = sim_out.p./(wages./prods);
        
        tc    = sum(costs.*sim_out.mu);
        num    = sum(costs.*sim_out.mu.*mkps);
        
        mu_cw = num./tc;
        
        figure(47)
        subplot(3,3,7)
        hold off
        plot((mu_cw)); title('Cost-weighted markup'); drawnow;
        xlim([1 T_plot]);

        figure(47)
        subplot(3,3,8)
        hold off
        plot((sim_out.W)); title('Wage');  drawnow;
        xlim([1 T_plot]);

%         % plot mass of entrants
         figure(47)
         subplot(3,3,9)
         hold off
         plot(sim_out.M);
         title('Mass of entrants');
         xlim([1 T_plot]);
         drawnow;
 
        diffC = sim_out.C - sim.C;
        diffC = abs(diffC)/max(abs(diffC));
        diffC = diffC;
        diffC = (1-options.damp_mit)*diffC;
        sim.C =(1-diffC).*sim.C + diffC.*sim_out.C;
        
        diffW = sim_out.W - sim.W;
        diffW = abs(diffW)/max(abs(diffW));
        diffW = diffW;
        diffW = (1-options.damp_mit)*diffW;
        sim.W =(1-diffW).*sim.W + diffW.*sim_out.W;
        
        
        cEupdate = .001;

         if options.shock == 'tfp'
             cEnew = param.ce.*sim.W;
             sim.cE = (1-cEupdate)*sim.cE + cEupdate*cEnew;
         end    
        
    end
    
    sim_out.entry = [];
    eq_mit = sim_out;
    
end